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ABSTRACT 


This  design  project  dealt  with  a  blending  problem  in  the  format  of  a  multivariable, 
constrained,  nonlinear  optimization  problem,  where  a  typical  application  of  this 
problem  may  be  an  incinerator  that  burns  waste.  The  incinerator  accepts  a  blend  of 
several  sources,  each  of  which  has  particular  characteristics  in  terms  of  thermal  output. 
Should  the  different  sources  be  inappropriately  blended,  it  could  result  in  a  safety 
hazard  for  the  entire  incinerator.  The  problem  is  to  determine  the  appropriate 
proportions  so  that  the  maximum  value  of  thermal  output  for  the  blend  is  minimized. 
The  maximum  value  will  be  approximated  by  the  mean  of  the  blend  plus  a  selected 
number  standard  deviations.  The  first  phase  of  the  project  was  a  brief  study  of  similar 
blending  problems  by  searching  relevant  literature.  It  was  found  that  very  little  has 
been  done  on  this  problem  in  terms  of  computer  implementation  for  industry. 

Secondly,  various  possible  methods  for  solving  the  blending  problem  were 
considered.  Among  these  methods  were  classical  calculus,  gradient  searches, 
exhaustive  searches,  and  iterative  direct  searching.  An  extension  of  an  iterative  direct 
search  is  a  Monte  Carlo  simulation,  where  a  statistical  approach  is  taken.  The  Monte 
Carlo  search  method  is  based  on  the  theory  that  the  optimization  problem  has  a 
distribution  of  answers  and  that  a  random  sample  of  those  answers  will  yield  an 
answer  in  the  lower  tail  to  within  a  certain  degree  of  "accuracy."  Then,  the  feasible 
region  of  answers  will  be  halved  and  another  random  sample  of  possible  answers  will 
be  taken.  The  search  proceeds  in  this  manner,  taking  random  samples  from  an  ever- 
decreasing  region  of  feasible  answers.  The  method  focuses  on  the  optimal  answer. 

The  third  aspect  of  the  project  was  implementing  the  Monte  Carlo  method  in 
computer  code.  The  program  was  written  in  Visual  Basic  4.0  for  Windows  in  order  to 
achieve  user-friendliness  and  a  graphical  user  interface.  The  program  contains  both  an 
optimization  fimction  and  a  "What-if"  analysis.  A  "User's  Manual"  was  also 
documented  to  assist  the  user. 
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M.I.S.E.  Design  Project  Proposal 

Essentially,  my  design  project  will  consist  of  a  blending  problem.  A  typical 
application  of  this  problem  to  the  industrial  domain  would  be  an  incinerator  that  burns 
waste.  The  incinerator  accepts  material  from  several  sources,  each  of  which  has 
particular  characteristics  in  terms  of  thermal  output.  The  thermal  output  of  each  source 
is  assumed  to  be  normally  distributed,  having  its  own  mean  and  standard  deviation. 
The  sources  will  be  blended  to  create  the  final  burning  material.  Should  the  different 
sources  be  inappropriately  blended,  it  could  result  in  a  safety  hazard  for  the  entire 
incinerator.  The  problem  is  to  determine  the  appropriate  proportions  so  that  the 
maximum  value  of  thermal  output  for  the  blend  is  minimized.  The  maximum  value 
will  be  approximated  by  the  mean  of  the  blend  added  to  a  selected  number  standard 
deviations.  For  example,  an  addition  of  three  standard  deviations  results  in  a  0.00135 
probability  that  the  maximum  value  will  be  exceeded.  As  the  number  of  sources 
increases,  the  problem  quickly  becomes  more  complex.  The  goal  of  the  project  is  to 
produce  a  working  model  computer  program  that  will  allow  user  interface  so  that  it 
will  calculate  the  appropriate  proportions  of  each  source  to  produce  an  optimally  safe 
blend  for  the  incinerator.  Additionally,  the  program  will  allow  the  user  to  play  "what- 
if"  games  in  order  to  test  different  blending  situations.  This  feature  will  be 
advantageous  when  the  blending  has  constraints. 

The  first  phase  of  my  project  will  be  a  brief  study  of  similar  blending  problems.  In 
talking  with  my  advisor.  Dr.  Hool,  it  was  discovered  that  very  little  has  been  done  on 
this  problem  in  terms  of  computer  implementation  for  industry.  This  phase  is  an 
attempt  to  determine  what  is  currently  available  in  this  field. 

Secondly,  a  mathematical  analysis  will  be  done  on  the  problem  to  determine  the  best 
solutions  to  the  problems,  as  the  number  of  sources  increases.  An  initial  hypothesis  is 
that  search  mediods  for  optimization  will  need  to  be  utilized  since  a  closed-form 
solution  is  impossible  and  an  exhaustive  computer  search  will  be  impractical. 
Additionally,  the  general  program  requirements  and  options  will  be  determined. 

The  third  aspect  of  the  project  will  be  implementing  the  working  solution  in 
computer  code.  The  initial  proposal  is  to  write  the  program  in  Visual  Basic  in  order  to 
achieve  user-friendliness  and  graphical  user  interface.  After  the  program  is 
flowcharted,  the  coding  will  be  done.  Naturally,  testing  and  de-bugging  will  also  be 
required.  Optionally,  a  short  "User's  Guide"  will  be  documented. 

Finally,  the  write-up  will  be  completed  in  order  to  ensure  thorough  documentation. 
As  of  now,  the  schedule  for  completion  date  is  June  10,  1997.  This  will  allow  almost 
three  months  for  the  review,  oral  examination,  processing  of  results,  and  payment  of 
the  graduation  fee  well  before  the  deadline  two-weeks  prior  to  graduation. 
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PROBLEM  ORIGIN 


This  problem  originated  from  a  real-world  problem  involving  hazardous  waste 
disposal  through  burning  of  hazardous  waste  material  in  an  incinerator.  Various 
gources  of  hazardous  waste  material  are  inputs  to  the  incinerator,  and  the  sources  are 
blended  into  a  mixture  prior  to  insertion  into  the  incinerator.  Each  source  contains  a 
distribution  of  Btu  (British  thermal  unit)  content  (i.e.  the  Btu  content  is  not  a  single 
value  for  a  source).  The  incinerator  cannot  be  safely  operated  if  the  maximum  Btu  of  a 
blended  input  exceeds  800  Btu  (the  "safe  limit").  The  problem  is  to  determine  how  to 
blend  the  sources  in  order  to  avoid  burning  material  that  has  a  maximum  Btu  value 
exceeding  the  "safe  limit"  of  the  incinerator.  Specifically,  the  problem  is  to  determine 
the  proportion  of  each  source  to  put  into  a  blend  so  that  there  is  a  small  probability  (e.g. 
0.001)  that  the  maximum  Btu  will  exceed  the  "safe  limit." 

The  blend  can  be  modeled  as  a  linear  combination  of  source  Btu's.  Each  source  has 
its  own  distribution  of  Btu  content.  The  mean  Btu  of  source  i,  i  =  1,  2,  ,  n,  is  |i, ,  and 

the  variance  of  Btu  content  of  source  i  is  a, 2.  If  the  proportion  of  source  i  that  is  blended 
into  the  final  mix  is  ki ,  then  the  Btu  content  of  the  blend  is  a  linear  combination,  B,  of 
the  sources.  If  Xi  is  the  Btu  content  of  source  i ,  where  Xi  is  a  random  variable,  then 

B  =  'Zk,X, 

i 

and  B  has  a  mean 

i 

and  a  variance 

i 

if  the  sources  are  independent.  The  independence  assumption  is  easily  met  in  the 
problem.  If  the  Btu  content  of  each  source  is  normally  distributed,  then  B  is  also 
normally  distributed.  If  not  all  of  the  sources  have  normally  distributed  Btu's,  then  B 
may  still  be  approximately  normally  distributed  if  enough  sources  are  blended  (as  a 
result  of  the  Central  Limit  Theorem). 
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If  B  can  be  assumed  to  be  at  least  approximately  normally  distributed,  then  for  any 
arbitrarily  selected  probability,  p,  where  p  will  be  small,  a  corresponding  "maximum" 
value  for  B,  defined  as  Mp ,  can  be  written  as 

where  Zp  =  standardized  normal  distribution  value  corresponding  to  the 
probability  p. 

Mp  is  the  (lOO)pth  percentile  of  B,  and  is  the  value  of  B  beyond  which  only  (100)p%  of 
the  blend's  Btu's  will  lie. 

If  the  "safe  limit"  of  the  incinerator  is  S,  then  the  blending  problem  becomes  one  of 
selecting  the  h  such  that 

M,<S 


subject  to  the  cor\straints 

XJt,=l 

1=1 

0  <ki<bi  Vi  =  !,•••,« 
where  h  =  1  are  default  values. 

In  particular,  the  "ultimate"  solution  of  the  blending  problem  requires  determination 
of  min(Mp)  =  Mp*.  If  Mp*  <  S,  then  the  blend  satisfies  the  incinerator's  "safe  limit."  Of 
course,  there  are  two  additional  situations  that  might  occur,  namely 


1.  Mp*  >  S,  in  which  case  a  blend  cannot  be  created  which  may  be  "safely" 
processed  by  the  incinerator. 

2.  There  are  multiple  Mp  (corresponding  to  various  ki  combinations)  for  which 

Mp<S. 

The  blending  problem  specific  to  diis  project,  and  any  blending  problem  of  similar 
character,  is  a  multivariable  constrained  nonlinear  optimization  problem.  The  decision 
variables  are  the  h,  and  the  term 
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creates  the  nonlmearity  of  the  objective  fimction  Mp.  Further,  several  variations  on  the 
problem  can  be  created  such  as  maximizing  a  minimum  value  of  B  such  that  the 
minimum  value  exceeds  a  lower  "safe  limit." 


PROBLEM  STATEMENT 

The  objective  of  this  study  was  to  select  an  appropriate  method  to  solve  the  blending 
problem,  and  then  produce  a  working-model  computer  program  to  execute  the  selected 
method  for  solving  the  problem.  The  computer  program  was  written  in  Visual  Basic  4.0 
for  Windows  and  performs  two  functions,  namely 


1.  Determines  the  optimal  h  (called  h*)  to  minimize  Mp.  The  program  can 
determine  if  the  blend  meets  the  "safe  limit"  {Mp*  <  S). 

2.  Allows  the  user  to  conduct  "what-if"  analyses  for  user-supplied  fc  values  in 
order  to  determine  if  the  resulting  blend  will  satisfy  Mp  <S. 

Because  the  number  of  sources  will  be  at  least  two  and  perhaps  as  many  as  ten,  the 
resulting  problem  is  the  following  multivariable  nonlinear  programming  problem: 


subject  to 

tk,-i 

1=1 

0  <ki<bi  Vi  =  n 


LITERATURE  SEARCH 

In  order  to  evaluate  what  is  currently  available  concerning  blending  problems,  a 
search  of  relevant  literature  was  conducted.  The  problem  being  investigated  is  a 
multivariable  constrained  nonlinear  optimization  problem.  Because  the  problem  is 
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nonlinear,  linear  programming  is  not  a  viable  option.  Most  nonlinear  optimization 
problems  concern  economics  or  inventory  problems.  The  most  preferred  method  of 
solving  nonlinear  optimization  problems  is  by  using  a  gradient  search  procedure, 
which  employs  the  partial  derivatives  of  the  objective  fvmction.  Optimization  methods 
will  be  more  fully  discussed  in  the  next  section. 

Nicholson  (1)  comments  that  in  the  industrial  world  a  gap  has  been  formed  between 
the  people  who  develop  optimization  methods  and  those  who  use  them.  Nicholson's 
premise  is  that  the  most  accurate  optimization  methods  will  be  useless  unless  they  are, 
if  not  theoretically  comprehended,  at  least  pragmatically  understood  by  the 
practitioner.  An  objective  of  this  project  was  to  implement  a  powerful  and  perhaps 
complex  optimization  method  and  make  it  useful  so  that  a  layman  can  easily  solve 
blending  problems. 

In  examining  Nicholson's  second  volume  of  applications,  along  with  other  industrial 
mixing  books,  it  was  foimd  that  most  mixing  problems  deal  with  complex  mixing 
formulas,  and  tend  toward  chemistry  problems.  Virtually  no  material  was  found 
relating  directly  to  solving  blending  problems. 

A  helpful  resource  that  was  discovered  was  Experiments  with  Mixtures  by  Cornell. 
However,  Cornell's  approach  toward  the  mixing  problem  adopts  design  of  experiments 
and  response  surface  metiiodology.  The  purpose  of  designing  a  complex  experiment  is 
to  construct  an  equation  that  adequately  describes  the  response  surface  (8).  Once  the 
response  surface  is  estimated,  it  is  then  optimized  by  selecting  the  appropriate  levels  of 
the  factors.  In  this  project,  the  response  surface  is  assumed  to  be  known  exactly  and  the 
real  problem  lies  in  selecting  the  levels  of  the  factors  (in  this  case,  h). 

POSSIBLE  SOLUTION  ALGORITHMS 

Multivariable  constrained  nonlinear  optimization  problems  can  be  solved  using 
various  methods.  Several  methods  that  were  considered  are  briefly  described  below. 
Note  that  use  of  a  linear  approximation  of  a  nonlinear  objective  function  is  less  accurate 
and  was  thus  not  considered. 
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1.  Analytical  (Calculus)  Method 

One  method  to  solve  a  nonlinear  optimization  problem  is  the  classical  calculus 
method.  This  involves,  in  a  multivariable  problem,  determining  the  first  partial 
derivative  of  the  objective  function  (Mp  in  this  case)  with  respect  to  each  variable  (i.e., 
fc).  Since  the  constraint  Zfc  =  1  is  used  in  the  blending  problem,  only  (n-1)  partial 
derivatives  must  be  determined.  Each  partial  derivative  must  be  set  equal  to  zero, 
resulting  in  a  system  of  equations,  and  the  h  then  solved  for  to  produce  a  global  or  local 
maximum,  minimum,  or  a  saddle  point  solution.  This  is  not  easy  in  the  general  case, 
and  in  particular  for  the  blending  problem.  For  the  blending  problem  each  partial 
derivative  is  a  nonlinear  function,  and  solutions  for  the  k  could  require  a  difficult 
search  process.  A  general  formula  for  the  (n-1)  first  partial  derivatives  is  shown  in 
Appendix  1,  and  an  example  for  n  =  2  sources  is  shown  in  Appendix  2.  Since  solution 
via  the  classical  calculus  method  is  complex,  this  method  can  seldom  be  used 
satisfactorily  on  even  xmconstrained  problems  (Nicholson  69).  In  addition,  determining 
whether  extreme  points  are  local  minima,  maxima,  or  saddle  points  can  be  difficult. 

2.  Gradient  Search 

When  the  system  of  equations  obtained  by  setting  the  partial  derivatives  equal  to 
zero  cannot  be  solved  analytically,  then  a  numerical  search  procedure  should  be  used 
(Hillier  514-5).  Using  the  objective  fimction,  Mp,  the  goal  is  to  reach  a  point  eventually 
where  all  of  its  partial  derivatives  are  (essentially)  zero  (Hillier  515).  The  gradient,  a 
vector  of  the  (n-1)  partial  derivatives,  can  be  calculated  and  used  to  determine  the 
direction  of  the  search.  However,  a  gradient  search  procedure  itself  must  employ 
another  search  procedure  to  determine  how  far  to  advance,  at  each  step,  in  the  direction 
of  the  gradient.  Therefore,  this  method  is  somewhat  complex  and  not  the  best  method 
to  use  when  attempting  to  optimize  with  a  computer.  An  example  of  this  is  shown  in 
Appendix  3. 
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3.  Failure  of  Calculus  -  Direct  Search 

This  type  of  search  method  differs  from  the  previous  method  because  direct  search 
procedures  work  directly  with  the  objective  fimction  and  avoid  the  difficulties  of 
calculating  the  gradients  of  functions  (Nicholson  84).  Advantages  of  this  type  of  search 
method  are  that  not  only  is  it  less  complicated  than  the  gradient  search,  but  it  also  can 
readily  cope  with  constraints  (Nicholson  84).  However,  a  disadvantage  is  that  with  an 
increasing  number  of  variables  (sources),  the  required  time  for  computation  may 
become  large  (Nicholson  84).  This  search  method  uses  a  base  point  and  perturbatioiis 
to  slowly  move  toward  the  optimal  point.  An  example  of  this  procedure  is  included  in 
Appendix  4. 


4.  Exhaustive  Search 

Since  the  previous  methods  are  somewhat  complex,  a  more  computer-adaptable 
method,  such  as  an  exhaustive  search,  was  considered.  This  method  employs  an 
examination  of  all  possible  combinations  of  fc  values,  and  then  identifies  the  set(s)  of  k, 
values  that  will  optimize  the  objective  function.  Depending  on  the  amount  of  accuracy 
desired  in  a  solution,  the  number  of  calculations  required  to  arrive  at  a  solution  can  be 
easily  determined.  For  example,  if  ±  1%  accuracy  is  desired  (i.e.,  h  =  0.00,  0.01,  ..., 
1.00),  then  each  source  (assuming  a  full  range  of  0.00  -  1.00)  would  have  101  different 
possible  ki  values.  The  following  table  shows  the  number  of  combinations  that  must  be 
examined  for  an  exhaustive  analysis  (assuming  n-1  independent  sources): 


Number  of  Sources 

Possible  Solutions 

2 

101 

3 

10, 201 

4 

1, 030, 301 

9 

1.08  X 10  16 

10 

1.09  X 10  18 

9 


However,  should  ±  0.01%  accuracy  be  desired,  each  source  will  have  10,001  possible 
values.  The  number  of  feasible  solutions  then  increases  considerably,  as  shown  in  the 
following  table  (assuming  n-1  independent  sources): 


Number  of  Sources 

Possible  Solutions 

2 

10, 001 

3 

100, 020, 001 

4 

1.00  X 10  12 

9 

1.00  X  10  32 

10 

1.00  X  10  36 

This  method  of  searching  has  the  disadvantage  of  being  computationally  intensive, 
even  for  a  fast  computer.  However,  there  are  several  advantages  of  this  method.  First, 
it  is  simple  to  execute.  A  computer  simply  uses  loops  to  search  over  the  entire  range  for 
each  variable,  comparing  the  current  solution  to  the  currently  optimal  solution. 
Second,  this  method  guarantees  optimality  within  a  certain  limit  of  accuracy.  Because 
every  possible  combination  is  considered,  optimality  can  be  assured.  It  is  the 
tremendous  speed  and  error-free  calculation  ability  of  the  computer  that  this 
exhaustive  search  method  seeks  to  exploit  (Conley  5). 

5.  Iterative  Search  (Strategic  Guessing) 

Although  the  exhaustive  search  will  work  in  all  cases  as  well  as  guarantee  optimality, 
it  must  be  admitted  that  much  of  the  computational  effort  amoimts  to  sheer  waste. 
Many  of  the  calculations  can  quickly  be  discarded  as  non-optimal.  However,  is  there  a 
way  for  the  computer  to  do  this?  Specifically,  if  an  algorithm  can  quickly  eliminate  a 
number  of  answers,  then  very  large  problems  can  be  solved  with  substantially  fewer 
calculations.  For  example,  as  indicated  earlier,  an  exhaustive  search  for  ten  variables 
ranging  from  0-100%  at  0.01%  increments  would  require  about  lO^^  calculations, 
disregarding  the  comparisons.  However,  what  if  an  algorithm  was  able  to  make  a  few 
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"smart"  decisions  and  reduce  the  number  of  calculations  to  ?  This  concept  of 
looking  only  at  strategic  points  and  using  those  values  to  make  decisions  is  the  driving 
force  behind  design  of  experiments  (DOE).  However,  DOE  can  be  used  to  adequately 
imderstand  the  response  surface,  then  assumes  that  optimization  will  easily  follow.  For 
this  project,  the  response  surface,  Mp,  is  known  and  its  optimization  is  the  problem. 
Nevertheless,  the  experimentation  methodology  may  be  used  to  strategically  locate  the 
optimal  point.  Cornell  explains  that  a  mixture  problem  with  n  ingredients  results  in  a 
factor  space  called  a  simplex.  The  geometric  description  of  the  factor  space  containing 
the  n  components  consists  of  all  points  on  or  inside  the  boundaries  (vertices,  edges, 
faces,  etc.)  of  a  regular  (n-l)-dimensional  simplex  (Cornell  6).  Three  sources  (n  =  3) 
result  in  a  simplex  factor  space  that  is  an  equilateral  triangle,  and  for  n  =  4  the  simplex 
is  a  tetrahedron  (Cornell  6).  Therefore,  the  next  step  is  to  determine  a  way  to 
selectively  examine  points  in  the  simplex;  one  possible  approach  is  to  evenly  divide  the 
simplex,  taking  observations  throughout.  An  ordered  arrangement  consisting  of  a 
uniformly  spaced  distribution  of  points  on  a  simplex  is  known  as  a  lattice  (Cornell  22). 
A  lattice  design  can  be  more  accurate  as  the  number  of  points  in  the  simplex  is 
increased.  A  lattice  design  is  said  to  have  degree  m  where  the  points  in  the  simplex  are 
spread  such  that  the  values  of  each  component  are 

1  2 

kj  =  0, — , — ,••',1 

m  m 

if  0  <  <  1.  Thus,  a  lattice  design  may  be  categorized  by  its  degree  and  by  the  number 

of  components  in  the  mixture.  The  number  of  selected  points  in  a  {n,m}  simplex-lattice 
design  will  be 

{n+m—  l)! 
m!(n-l)!  ' 

For  example,  a  {3,2}  design  will  have  6  points,  while  a  {10,4}  design  will  have  715 
points.  A  table  showing  the  total  number  of  points  for  a  [n,m]  lattice  is  shown  in 
Appendix  5.  Then,  a  further  lattice  could  be  used  on  a  ,  selected  region  of  the  first 
simplex,  so  that  the  optimal  point  is  slowly  brought  to  focus.  An  example  with  n  =  2  is 
shown  in  Appendix  6. 
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6.  Iterative  Search  (Monte  Carlo  Simulation) 

This  was  the  method  selected  for  use  with  the  blending  problem.  Instead  of 
attempting  to  equally  divide  the  simplex  into  regions  as  is  done  with  an  iterative  search 
method,  a  Monte  Carlo  simulation  method  examines  a  random  sample  of  points  inside 
the  simplex  (Conley  20).  Based  on  the  optimal  value  in  the  sample,  the  simplex  region 
can  be  halved  and  another  random  sample  of  points  can  be  taken.  The  computer  can 
select  the  optimal  sample  value  at  each  iteration  to  be  the  centroid  for  the  next 
iteration's  simplex,  halving  the  simplex  region  as  well.  This  technique  employs 
statistical  probability  and  then  relies  on  the  massive  speed  of  the  computer  to 
determine  an  optimal  solution. 

As  with  any  optimization  technique,  both  advantages  and  disadvantages  exist.  A 
primary  advantage  is  that  the  method  uses  far  fewer  calculations  than  an  exhaustive 
search.  Conley  stated  (as  of  1981)  that  a  computer,  using  Monte  Carlo  techniques, 
found  a  true  optimum  solution  to  a  problem  with  lO^o  feasible  solutions  in  about  a 
minute  on  a  medium-sized  computer,  where  an  exhaustive  search  technique  on  a 
problem  of  this  scale  would  take  the  computer  about  32  million  trillion  years  (Conley 
20).  Other  advantages  include  the  ability  to  use  an  accurate  (nonlinear)  model  with  no 
concern  for  linearity,  continuity,  or  differentiability  assumptions.  One  disadvantage  is 
that  optimality  cannot  be  assured  because  the  method  is  based  on  randoirmess,  and  the 
method  could  conceivably  miss  an  optimal  point.  However,  the  probability  that  a 
Monte  Carlo  solution  investigating  1,000,000  combinations  will  come  within  0.00001  of 
the  optimal  solution  is  1  -  0.99999^'°®°-°°°  =  0.9999546.  Therefore,  the  probability  is 
almost  overwhelming  that  at  least  a  near  optimum  solution  will  be  foxmd. 

The  basis  for  this  method  is  that  all  feasible  solutions  map  a  distribution  of  objective 
function  values  (see  Appendix  7).  Most  values  tend  toward  the  objective  function 
mean,  and  the  optimal  solution  (to  the  specific  problem  of  this  project,  Mp)  is  at  the  left 
tail  of  the  distribution.  At  each  iteration  of  the  simulation,  a  specified  number  of 
samples,  say  4000,  are  selected.  The  probability  that  one  random  sample  will  not 
produce  a  value  in  the  lower  0.01  tail  of  the  Mp  distribution  is  0.99.  However,  the 
probability  that  none  of  the  4000  samples  will  produce  a  value  in  the  lower  0.01  tail  of 
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the  distribution  is  0.99  ^ooo  ^  which  is  about  3.5  x  10  -is.  With  each  iteration,  the  feasible 
region  is  made  smaller  and  smaller,  so  that  another  4000  samples  will  most  likely  yield 
a  solution  in  an  even  lower  tail  of  the  distribution,  say  0.0001.  After  several  iterations 
(<  10),  the  optimal  solution  is  virtually  guaranteed.  Therefore  the  optimal  solution  will 
most  likely  be  found  after  the  computer  evaluates  a  mere  4000  x  10  =  40,000  possible 
solutions,  which  is  a  drastic  reduction  in  calculations  compared  to  the  exhaustive 
search  method.  The  only  situation  that  would  disrupt  the  Monte  Carlo  approach 
would  be  a  distribution  (in  an  optimization  problem)  that  had  one  isolated  optimal 
solution  that  was  separate  from  the  remainder  of  the  solutions  (Conley  189,  see 
Appendix  8). 

The  optimization  procedure  works  as  follows: 

1.  Assign  the  initial  optimal  solution  by  letting  fc  =  1/n  for  z  =  1,  2,  . . .  ,  n,  where 
Iki  =  1  and  0  <  h  <  h  .  Call  this  solution  h*.  Its  corresponding  objective 
function  value,  Mp,  will  be  Mp*. 

2.  Set  iteration  number,  /  =  1- 

A.  Randomly  generate  4000  samples,  where  for  sample  s,  s  =  1,  2,  ...  , 
4000,  the  ki  are  randomly  selected  on  the  interval  0  <ki<  h  and  where 
Iki  =  1.  Compute  the  corresponding  Mp  value  for  each  sample. 

B.  For  each  sample  s  =  1,  2,  . . .  ,  4000,  compare  Mp  with  Mp*.  If  Mp  <  Mp*, 
then  replace  Mp*  with  Mp  and  let  h*  =  fc. 

C.  Using  ki*,  z  =  1, 2, ... ,  n,  as  the  centroid  for  source  z,  determine  the  new 
reduced  range  of  each  h  for  the  next  iteration  where 

max  {ki*  -  2-i-b  0)  <fc  <min  (fc*  +  2-i-b  h). 

3.  For  iterations  j  =  2,3, ... ,  10,  use  h*  and  Mp*  from  the  previous  iteration. 

A.  Randomly  generate  4000  samples,  where  for  sample  s,  s  =  1,  2,  ...  , 
4000,  the  ki  are  randomly  selected  on  the  interval 

max  {ki*  -  2-i- 1, 0)  <  ki  <  min  {h*  +  2-i  •  b  bi) 
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obtained  from  the  previous  iteration,  where  Sfc  =  1  and  a 
corresponding  Mp  is  computed. 

B.  For  sample  s  =  1,  2,  ,  4000,  compare  Mp  with  Mp*.  If  Mp  <  Mp*,  then 

replace  Mp*  with  Mp  and  let  fc*  =  h. 

C.  Using  ki*,  i  =  1, 2, ... ,  n,  as  the  centroid  for  source  i,  determine  the  new 
reduced  range  of  each  k  for  the  next  iteration  where 

max  (fc*  -  2-i-i,  0)  <  k  <  min  {k*  +  2-i-\  bi) 
where  j  indicates  the  iteration  just  completed. 

D.  Return  to  step  3  imtil  10  iterations  have  been  completed. 

4.  The  identified  optimal  solution  is  k*,  with  an  objective  function  value  of  Mp*. 

For  example,  assume  that  the  problem  consists  of  three  sources,  that  the  optimal 
solution  is  0.1,  0.3,  and  0.6  for  each  source  respectively,  resulting  in  art  objective 
function  value  of  50,  and  that  each  source  has  boimds  placed  at  0.0  and  0.8: 

Step  1.  The  initial  solution  (fc*)  would  be  set  at  0.3333  for  each  source,  with  an 
objective  function  value  of,  say,  Mp*  =  62. 

Step  2  (A)  The  algorithm  would  generate  the  first  sample  of  proportion  values, 
say  ki  =  0.2289,  ki  =  0.5372,  and  ks  =  0.2339,  and  (B)  evaluate  the  objective  function 
to  be  Mp  =  63.87.  Therefore,  the  second  random  sample  would  be  generated  and 
again  compared  to  the  initial  optimal  solution.  Assume  that  after  4000  samples, 
the  best  solution  thus  far  is  fa*  =  0.1019,  fa*  =  0.2901,  and  fa*  =  0.6080,  with  an 
objective  fimction  value  of  Mp*  =  52.17.  (C)  This  solution  would  form  the  basis 
for  determining  the  new  range  of  possible  values.  For  iteration  1,  the  algorithm 
attempts  to  add/subtract  2-2  =  0.25  from  each  fa*  value  to  obtain  the  new  reduced 
range.  The  first  source  would  now  have  a  range  of  (max  (0.1019  -  0.25,  0.00))  = 
0.00  and  (min  (0.1019  +  0.25,  0.8000))  =  0.3519,  the  second  a  range  of  0.0401  - 
0.5401,  and  the  third  a  range  between  0.3580  -  0.8000.  Thus,  the  second  source 
has  a  full  range  of  0.5000  while  the  other  two  were  limited  by  their  boimds, 
resulting  in  smaller  ranges. 
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Step  3  (A/B)  Iteration  2  would  then  begin  and  another  4000  samples  would  be 
generated  within  this  smaller  range  of  ki  values.  (C)  Assume  this  time  that  after 
the  4000  additional  samples,  the  best  proportions  found  thus  far  are  0.1009, 
0.2987,  and  0.6004,  with  an  objective  fimction  value  of  Mp*  =  50.19.  This  solution 
would  form  the  basis  for  determining  the  new  range  of  possible  values.  For  each 
source,  the  algorithm  now  attempts  to  add/ subtract  2  ■  ^  =  0.125  from  each  h* 
value  to  obtain  the  new  reduced  range.  The  first  source  would  now  have  a  range 
of  0.0000  -  0.2269,  the  second  a  range  of  0.1737  -  0.4237,  and  the  third  a  range  of 
0.4754  -  0.7254.  This  time,  the  second  and  third  sources  have  a  maximum  range 
of  0.25  while  the  first  has  a  smaller  range.  (D)  The  second  iteration  would  be 
complete  and  eight  more  iterations  would  follow,  each  taking  4000  samples  from 
an  ever-decreasing  range.  At  iteration  10,  the  range  for  each  of  the  sources  would 
be  2  -  9  =  1.91  X  10  -  3.  This  means  that  4000  random  samples  would  be  selected 
from  this  range  on  each  of  the  h  values,  which  would  guarantee  optimality  with 
acceptable  accuracy. 

Step  4  The  optimal  solution  of  h*  =  (0.1000,  0.3000,  0.6000)  for  i  =  1,  2,  3,  with  an 
objective  function  value  of  Mp*  =  50.00000  is  reported  as  the  final  answer. 

At  each  iteration,  the  feasible  region  is  at  least  halved  and  most  likely  reduced  much 
more.  If  j  is  the  iteration  number,  where  ;  =  1,  2,  ...,  10,  and ;  =  1  (first  iteration),  with 
ki*  =  0.7  and  the  bounds  0.00  <  ki  <  0.90,  then  for  the  next  iteration  (i.e.,  j  =  2), 
max  (0.7  -  2  - 1  - 1, 0.00)  <  ki  <  min  (0.70  +  2  ■  i  -  b  0.90) 
max  (0.7  -  0.25, 0.00)  <  ki  <  min  (0.70  +  0.25, 0.90) 
max  (0.45, 0.00)  <  ki  <  min  (0.95, 0.90) 

0.45  <ki<  0.90. 

Essentially,  the  procedure  uses  (n-l)-dimensional  "rectangles"  that  search  through 
the  n  dimensional  space  (always  staying  inside  the  constraints)  imtil  the  rectangles 
focus  on  the  optimal  answer  (Conley  229).  According  to  the  laws  of  probability,  the 
algorithm  will  quickly  find  the  optimal  solution. 
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The  Monte  Carlo  method  disregards  classical  optimization  methods  in  favor  of 
relying  on  the  pure  speed  of  the  computer  combined  with  the  laws  of  probability.  For 
any  optimization  problem  of  one  himdred  variables  or  less,  Conley  (245)  recommends 
using  this  method  to  optimize  a  function. 

MONTE  CARLO  SEARCH  COMPUTER  PROGRAM 
The  program  has  two  options:  1)  optimize  a  blend,  and  2)  user  "What-if"  analysis. 

The  first  option  will  allow  the  user  to  input  the  following  values: 

1.  Number  of  sources,  n  (selected  from  the  range  2-10) 

2.  For  each  source: 

A.  Source  mean,  |Xi 

B.  Source  variation :  variance,  or  standard  deviation,  <Ji 

C.  Source  upper  boimd,  h  (1.00  =  default  value);  (bn  must  remain  at  1.00) 

3.  "Safe  limit,"  S  (used  for  comparison) 

4.  Probability  p  (proportion  of  actual  blend  output  that  will  exceed  "maximum" 
value,  Mp*.  The  program  will  then  determine  the  appropriate  number  of 
standard  deviations,  Zp,  to  add  to  the  expected  mean  of  the  blend  to 
correspond  with  this  probability.  Note  that  possible  p  values  are:  0.25,  0.10, 
0.05, 0.025, 0.01, 0.005, 0.0025, 0.001, 0.0001). 

The  program  will  output: 

1.  The  expected  output  of  the  blend,  pe 

2.  The  standard  deviation  of  the  blend,  Ob 

3.  The  optimal  minimum  "max"  value,  Mp* 

4.  A  comparison  of  optimal  and  "Safe  limit"  (SAFE /UNSAFE) 

5.  The  corresponding  ki*  for  each  source 

The  second  option  of  the  program  will  be  for  user  "What-if"  analysis.  The  user 
inputs  will  be: 

1.  Number  of  sources,  n  (selected  from  the  range  2-10) 

2.  For  each  source: 
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A.  Source  mean,  [li 

B.  Source  variation:  variance,  or  standard  deviation,  Oi 

C.  Proportion  of  the  source,  h,  in  the  blend  (in  keeping  with  Zfc  =  1) 

3.  "Safe  limit"  value,  S 

The  program  will  then  give: 

1.  The  expected  output  value  of  the  given  blend,  |Xb 

2.  The  standard  deviation  of  the  blend,  Ob 

3.  Proportion  of  the  blend  output,  B,  that  will  exceed  the  "Safe  limit." 

VISUAL  BASIC  4.0 

Microsoft's  Visual  Basic  4.0  for  Windows  was  the  programming  language  selected  to 
code  the  computer  program.  "Visual  Basic  is  an  object-oriented /event-driven 
programming  language  that  is  easy  enough  for  a  nonprogrammer  to  use,  yet 
sophisticated  enough  to  be  used  by  professional  programmers"  (Zak  5).  Although 
Visual  Basic  is  a  relatively  new  language,  several  advantages  that  the  language  offers 
made  this  selection  practical.  Among  the  advantages  are: 

•  Visual  Basic  programs  are  rim  in  a  Windows  environment,  making  the  user 
interface  simple,  neat,  and  user-friendly 

•  Visual  Basic  programs  can  be  stored  in  an  executable  file  (*.exe)  so  that  any 
user  with  a  Windows  95  operating  system  may  run  the  program 

•  Visual  Basic  allows  the  programmer  to  spend  more  time  in  coding  the 
program  details  while  user  interface  is  virtually  self-coded  by  the  program 
itself. 

A  Visual  Basic  program  is  structured  around  forms  and  objects.  Forms  are  the 
different  windows  that  will  allow  user  interaction,  and  the  objects  are  the  various  text 
boxes,  drop  menus,  and  click  commands  that  appear  on  each  form.  Therefore  the  first 
step  in  Visual  Basic  programming  is  to  develop  the  user  interface  in  the  framework  of 
forms  and  objects.  After  this  major  step,  instructions  concerning  each  object  on  the 
form  are  given  as  code. 
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This  program  has  four  basic  forms.  The  first  form  is  the  "Control  Panel."  From  this 
form,  the  user  may  move  to  one  of  the  major  program  functions  or  exit  the  program. 
The  second  major  form  is  the  "Optimization"  form.  This  form  contains  the  overall 
inputs  and  all  outputs  for  the  optimization  problem.  The  third  form  is  the  "What-if" 
analysis  form  which  contains  overall  inputs  and  all  outputs  for  the  "What-if"  analysis. 
The  last  form  actually  contains  a  sub-form  for  each  source.  These  sub-forms  allow  the 
user  to  input  data  specific  to  each  blending  source. 

In  order  to  allow  the  reader  and  perhaps  others  who  wish  to  expand  on  this  project  to 
more  fully  imderstand  the  program  as  a  whole  and  the  optimization  procedure, 
pseudocode  is  given  as  Appendix  9.  Furthermore,  the  actual  code  for  the  optimization 
of  the  blending  problem  with  two  sources  is  given  as  Appendix  10,  while  Appendix  11 
contains  the  code  used  in  the  "What-if"  analysis  for  two  sources.  The  reader  will  note 
that  all  Visual  Basic  4.0  "comments"  are  preceded  by  an  apostrophe  (')• 

PROGRAM  VALIDATION 

In  order  to  verify  whether  Blend  was  properly  optimizing  proportions,  nine  sample 
problems  (one  for  each  n  =  2,.. .,10)  were  optimized.  Before  optimizing  with  Blend,  the 
problems  were  solved  with  LINGO,  a  nonlinear  optimization  software  package.  Each 
problem  was  rim  five  times  while  comparing  only  the  Mp*  values,  since  it  is  true  that 
the  actual  proportions  can  vary  substantially  while  maintaining  extremely  similar  Mp* 
values.  From  all  45  runs,  the  average  deviation  from  the  optimal  answer  was  0.66%. 
Of  interest  was  the  fact  that  most  of  the  error  came  from  one  problem  that  Blend 
seemed  to  have  more  difficulty  in  optimizing,  while  no  unique  features  of  the  problem 
could  be  established.  Overall,  Blend  demonstrated  a  strong  optimization  potential; 
experimental  results  for  problems  where  n  =  3, 7,  and  10  are  shown  in  Appendix  12. 
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APPENDIX  1 

General  Formula  Using  (w-1)  Partial  Derivatives 


General  Formula  Using  (n-1)  Partial  Derivatives 


For  this  derivation,  assume: 


Source 

Proportion 

Mean 

Standard  Deviation 

1 

ki 

SiJUi 

riai 

2 

ki 

SiMi 

rzCTi 

n-l 

kn-l 

Sn-lJUl 

rn-l<Jl 

n-l 

n 

1-Z^, 

i=l 

SnMl 

TnCTl 

and  where  Si  =  ri  =  1.  Let  the  subscript  B  denote  the  blend. 
Blend  Mean: 
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Blend  Standard  Deviation: 
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Now,  define  Mp  =  //b  +  Z<tb 


n-l 

'n-l  f  ^ 

i" 

f^P=Mi 

+  Z(j^ 

i=\  \  i=i  ^ 
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Taking  the  partial  derivative  results  in  (n-1)  equations: 


n-\  f  «-l  'N  ^ 

,=1  V  i=i  J 


-K 


n-1 


Ikjvj  +lk^rl-lrl  +2r^\ 


y/  =  1,2,. 


y  V 

(5  -5j  +  Z<T, 

Vi  ^''n 

\  J  "/  ‘ 

V  V 

n-1 


i=I  > 

Vy  =  l,2,...,n-1 


n-1 


+  rj 


1=1 


1-Z^i 

i=i  / 


-ly 


r  z^"-* 

n-1 

y  n-1 

2" 

-O+^O-I 

V/+''»|Z^--M 

1-Z*,, 

^  v.,=i  y; 

1=1 

\  ,=1  / 

_ 

,-iy 


y/  =  i,2,...,«-i 


^p 


=  A  (^y  -  -^n)  +  ^^1  V/  -  '■«  1 1  -  Z 

V  ^  i=i 


n-l  W 


JJ 


n-\  f  n-\  rt-1  n-l 

Z^iV  +'•»  I  i-2Z^,-  +zz^,^ 

^  1=1  1=1  m=l  ^ 


i=l 


y/  =  l,2,...,n-l 
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APPENDIX  2 


Example  of  an  Analytical  Solution  when  n  —  2  Sources 


Example  of  an  Analytical  Solution  when  n  =  2  Sources 


Given : 

Hi  =  50  Oi  =  2  Hz  =  40  02  =  6 

assume 

Za  =  0.00135,  so  Z  =  3.00 

then 

Si  =  l  ri=l  S2=.8  r2  =  3 


solution  will  require  (n  -1)  =  1  derivative  : 


n-\ 

,=1 


n-1 


^  n— 1  /I— 1  n-~\ 


H  1  -  2E  +  Z  Z 

V  i^\  ;„=t 


2=1 


r'A 


=  50  (l  -0.8  )  +  3-2  (A:,-f-3'(l-)t,))[A:,-l"+3'(l-2^,+A:f)] 


set  derivative  equal  to  zero  and  solve  for  k\ : 

m 


=  10  +  (60A:,  -54)[l0yt/  -  ISAj  +9]  ^  =  0 


-  lOyllOk^  -m,+9  =  60k,  - 54 


100  ( 10fci2  -  ISfci  +  9)  =  3600fci2  -  6480fc  +  2916 
26fci2-46.8fci  + 20.16  =  0 


using  the  quadratic  formula ; 


k,= 


46.8  ±  V46.8^-4-26-2016 


2-52 


fci  =  1.086  or  0.7139 


then 


fci  =  0.7139  and  kz  =  0.2861 
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APPENDIX  3 


Example  of  Gradient  Search  Procedure  {n  =  2) 


Example  of  Gradient  Search  Procedure  ( w  =  2 ) 

(adapted  from  Hillier/Lieberman) 


Given : 

|ii  =  50 
assume 

8  =  .005 

then 

Si=  1 


oi  =  2  |j,2  =  40  02  =  6 

Za  =  0.00135,  so  Z  =  3.00 

ri=l  S2=.8  r2  =  3 

Mp  =  lOit,  +40+67l0it,^-18A:,+9 

=  10  +  (60A,  -  54X1  Oit,'  - 1 8A:,  +  9)'^ 


begin  at  k\  =  0.5 


dMp 

dk^ 


•  =  -0.7331 


VM^  =  (-0.7331) 

fci=  0.5  -  f  (-0.7331)  =  0.5  +  0.7331f 


M.  =  1 0(0.5  +  0.733 10  +  40  +  6^10(0.5  +  0.73310"  - 18(0.5+  0.733 h)  +  9 


M  =  45  +  7.33 \t  +  6V2.5  +  7.33 1/  +  5.3744/^  -  9  - 13.1958/ +9 


Mp=AS  +  733  It  +  6V5.3744^"  -5.8648^  +  2.5 


one-dimensional  search  to  find  min  (f  >0)  t*  of  Mp 


dMp 

dt 


=  7.331  + 


6- 0.5(10.7488/ -5.8648) 
75.3744/^-5.8648/ +2.5 


=  0 


32.2464/- 17.5944 


=  -7.331 


75.3744/^ -5.8648/ +  2.5 
(322464/-17.5944)"  =(-7.33l)"(5.3744/^  -5.8648/  +  2.5) 
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1039.8303/^ -1134.7121/  +  309.5629  =  288.8396/^-315.1955^  +  134.359 


750.9907/^  -819.5166/  + 175.2039  =  0 
t  =  0.7994  or  t  =  0.2918 


U  =  0.5  +  0.7331(  0.7994 )  or  fc  i  =  0.5  +  0.7331(  0.2918 ) 
=  1.086  =  0.7139 

*  not  feasible  *  feasible 


<M  11166 

(0.7 139)  =  ^  =10-  -j—  =  -0.002 

^  dk,  Vl.2463 

which  is  sufficiently  close  since  e  =  .005. 
Therefore,  ki  =  0.7139  and  fo  =  0.2861. 
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Example  of  Direct  Search  Procedtire 


Example  of  Direct  Search  Procedure 

(adapted  from  Nicholson) 


Given : 

[4,1=50  ai=2  [42  =  40  02=6 

assume 

s  =  .01  Z  =  3.00 

then 

Si=l  ri=l  S2=  .8  r2  =  3 


Mp  =  lOA:,  +40+6yll0k^ -lSk,+9 


begin  at  fci  =  0.5 

perturbation  step  size,  6  =  0.01 

*Base  point,  =  0.5 

Mp  ( BW )  =  54.4868 

Mp  ( BW  +  6)  =  54.4357  accept  and  let 

ToO)  =  0.51 

*  check  "double"  step 

Mp  (B(o)  +  28)  =  54.3861  accept 

TiO)  =  0.52 

*  check  "double"  step 

Mp  (B(“)  +  2(25))  =  54.2913  accept 

T2(i)  =  0.54 

*  check  "double"  step 

Mp  (B(o)  +  2(45))  =  54.1225  accept 

T3(i)  =  0.58 

*  check  "double"  step 

Mp  (BW  +  2(85))  =  53.8894  accept 

T4(2)  =  0.66 

*  check  "double"  step 

Mp  (B(o)  +  2(165))  =  54.0910  reject 

so  T4(«  =  0.66  =  B(i) 

*  begin  again 

Mp  (BW  +  5)  =  53.8724  accept 

To(2)  =  0.67 

*  check  "double"  step 

Mp  (B(^)  +  25)  =  53.8586  accept 

Ti(2)  =  0.68 

*  check  "double"  step 

Mp  (B(i)  +  2(25))  =  53.8411  accept 

It 

o 

o 

*  check  "double"  step 

Mp  (B(i)  +  2(45))  =  53.851  reject 

so  T2(2)  =  0.70  =  B(2) 

,'y 
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*  begin  again 

*  check  "double"  step 


Mp  (B(2)  +  5)  =  53.8377  accept 
Mp  (B(2)  +  25)  =  53.8381  reject 


so  ToP)  =  0.71  =  BP) 


ToP)  =  0.71 


* 


* 


begin  again 


Mp  (BP)  +  5)  =  53.8381  reject 

since  no  step,  try  step  back 

Mp  (BP)  -  5)  =  53.8411  reject 

since  step  in  either  direction  is  worse,  take  current  base  as  optimal 
B(3)  =  0.71  so  h  =  0.71  =  0.29 
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Total  Points  of  a  Lattice  Design 


Total  Points  of  a  Lattice  Design 
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Example  of  a  Lattice  Design  Search  Procedure  ( «  -  2 ) 


Example  of  a  Lattice  Design  Search  Procedure  ( w  =  2 ) 


Given : 

|j,i  =  50  Oi  =  2  |j,2  =  40  02  =  6 

assume 

8  =  .005  Z  =  3.00  lattice  order  =  5 

then 

si=l  ri=l  S2=.8  r2  =  3 


Mp  =  lOA:,  +40+6yll0k^-m,+9 


*  search  for  optimal  h  begiiming  with  endpoints  at  (0,1) 


Mp(0)  = 

58 

Mp(0.25)  = 

56.08 

Mp(0.50)  = 

54.49 

Mp(0.75)  = 

53.86 

Mp(l)  = 

56 

*  re-center  search  and  continue  . . . 


Mp(0.50)  = 

54.49 

Mp(0.625)  = 

53.97 

Mp(0.75)  = 

53.86 

Mp(0.875)  = 

54.46 

Mp(l)  = 

56 

*  re-center  search  and  continue  ... 


Mp(0.625)  = 

53.97 

Mp(0.6875)  = 

53.850  * 

Mp(0.75)  = 

53.86 

Mp(0.8125)  = 

54.054 

Mp(0.875)  = 

54.46 

*  re-center  search 

and  continue 

Mp(0.625)  = 

53.97 

Mp(0.65625)  = 

53.8966 

Mp(0.6875)  = 

53.850 

Mp(0.71875)  = 

53.8378  * 

Mp(0.75)  = 

53.86 
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*  re-center  search  and  continue  . . . 


Mp(0.6875)  = 

53.850 

Mp(0.703125)  = 

53.83959 

Mp(0.71875)  = 

53.8378  * 

Mp(0.734375)  = 

53.8457 

Mp(0.7500)  = 

53.86 

*  re-center  search  and  continue 

Mp(0.703125)  = 

53.83959 

Mp(0.7109375)  = 

53.83753  * 

Mp(0.71875)  = 

53.8378 

Mp(0.7265625)  = 

53.8405 

Mp(0.734375)  = 

53.8457 

*  re-center  search  and  continue 

Mp(0.703125)  = 

53.83959 

Mp(0.70703125)  = 

53.83827 

Mp(0.7109375)  = 

53.83753 

Mp(0.71484375)  = 

53.8373702 

Mp(0.71875)  = 

53.8378 

*  re-center  search  and  continue 

Mp(0.7109375)  = 

53.83753 

Mp(0.712890625)  = 

=  53.8373763 

Mp(0.71484375)  = 

53.8373702 

Mp(0.716796875)  = 

=  53.8375129 

Mp(0.71875)  = 

53.8378 

accuracy  currently  at  0.001953125,  which  is  sufficient  since  s  — 

fci  =  0.71484375  and  ki  =  0.28515625 


.005. 
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APPENDIX  7 

Distribution  of  Objective  Function  Values 


Distribution  of  Objective  Function  Values 

(taken  from  Conley) 


/i 
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APPENDIX  8 

Disruption  of  a  Monte  Carlo  Approach 


Disraption  of  a  Monte  Carlo  Approach 


one  extra, 
isolated  value 
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Pseudocode 


PSEUDOCODE 


This  appendix  gives  the  general  code  for  each  object  of  each  form  (i.e.  how  the 
program  will  respond  to  any  action  that  the  user  might  take).  Note  that  dimming 
involves  leaving  a  label  or  text  box  in  sight  while  not  allowing  the  user  to  access  it 
Also,  the  mouse  "focus"  is  defined  as  the  current  cursor  position  on  a  form. 

1.  Control  Panel 

When  the  user  clicks  on  "Optimize" 

Remove  the  "What-if '  form  from  view,  should  it  currently  be  active 
Show  the  "Optimize"  form  while  also  showing  the  Source  Data  form  minimized 
When  the  user  clicks  on  "What-if' 

Remove  the  "Optimize"  form  from  view,  should  it  currently  be  active 
Show  the  "What-if'  form  while  also  showing  the  Source  Data  form  minimized 
When  the  user  clicks  on  "Exit" 

Exit  the  program 


2.  Optimize 

When  the  form  loads 

Initialize  the  "Already"  boolean  variable  to  false 
Load  the  choices  into  the  "Number  of  Sources"  list  box 
Initialize  the  box  set  to  Number  Sources  =  2 
Load  the  choices  into  the  "Probability"  list  box 
Initialize  the  box  set  to  Probability  =  0.01 
Call  the  "Dim  Blend  Proportions"  procedure 
Call  the  "Un-dim  Bounds"  procedure 
When  the  user  changes  the  Number  of  Sources  using  the  list  box 
Hide/ Disable  Source  Data  forms  and  Optimal  Proportions  (sources  3-10) 
Un-dim  all  upper  bounds  for  sources  3-9 
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Show  appropriate  Source  Data  forms  and  Optimal  Proportions  (sources  3-10) 
Show  warning  message  about  forms  8-9-10  if  necessary  (Already  is  false)  and  set 
Already  to  true 

When  the  user  clicks  on  the  "Cancel"  button 
Close  Source  Data/  Optimize  forms 
Show  the  Control  Panel  form 
When  the  user  clicks  on  the  "Optimize"  button 
Declare  local  variables 
Initialize  variables 

Call  appropriate  optimization  procedure  (depending  on  number  of  sources) 
Display  output  values 
When  the  user  changes  the  "Safe  limit" 

Error-check  to  ensure  that  a  number  was  input 

Show  warning  message  if  user  did  not  input  a  numerical  value 

Return  mouse  focus  to  the  "Safe  limit"  text  box 

3.  "What-if'  Analysis 

When  the  form  loads 

Initialize  the  "Already"  boolean  variable  to  false 
Load  the  choices  into  the  "Number  of  Sources"  list  box 
Initialize  the  box  set  to  Number  Sources  =  2 
Call  the  "Un-dim  Blend  Proportions"  procedure 
Call  the  "Dim  Bounds"  procedure 
When  the  user  changes  the  Number  of  Sources  using  the  list  box 
Hide/ Disable  Source  Data  forms  (for  sources  3-10) 

Show  appropriate  Source  Data  forms  (for  sources  3-10) 

Show  warning  message  about  forms  8-9-10  if  necessary  (Already  is  false)  and  set 
Already  to  true 
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When  the  user  clicks  on  the  "Cancel"  button 
Qose  Source  Data  /  What-if  Analysis  forms 
Show  the  Control  Panel  form 
When  the  user  clicks  on  the  "Compute"  button 
Declare  local  variables 
Initialize  variables 

Call  appropriate  compute  procedure  (depending  on  number  of  sources) 

Display  output  values 

When  the  user  changes/finishes  changing  the  Maximum  Safety  Value  text  box 
Error-check  the  input  to  ensure  that  it  is  a  numerical  value 
Return  the  mouse  focus  back  to  the  Max  Safety  Value  text  box 

4.  Source  Data  (10  forms,  aU  identical) 

When  the  form  loses  the  mouse  focus 
Declare  local  variables 

Error-check  all  source  data  (that  all  values  are  numerical) 

Display  an  appropriate  message  if  error  exists 
When  the  user  changes/finishes  changing  any  of  the  following:  blend  proportion,  mean  output, 
output  variation,  upper  bound,  lower  bound 
Declare  local  variables 

Error-check  the  appropriate  data  (to  ensure  that  it  is  a  numerical  value) 

Check  to  ensure  that  numbers  are  feasible 

Display  error  message  if  necessary 

Return  mouse  focus  to  text  box  found  in  violation 

5.  Procedures 

Dim  Blend  Proportions  Procedure 

Dim  the  blend  proportion  label  and  text  box  on  each  Source  Data  form 
Dim  Bounds  Procedure 

Dim  the  bound  frame,  labels,  and  text  boxes  on  each  Source  Data  form 
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Un-dim  Blend  Projfortions  Procedure 

Un-dim  the  blend  proportion  label  and  text  box  on  each  Source  Data  form 
Un-dim  bounds  Procedure 

Un-dim  the  bound  frame,  labels,  and  text  boxes  on  each  Source  Data  form 
Max  Procedure 

Set  the  maximum  value  to  the  first  element  in  the  array 

Compare  each  element  in  the  array,  setting  it  to  maximum  value  if  it  is  larger 
Min  Procedure 

Set  the  minimum  value  to  the  first  element  in  the  array 

Compare  each  element  in  the  array,  setting  it  to  minimum  value  if  it  is  smaller 
Z- Alpha  Procedure 

Select  the  appropriate  case  of  the  user-defined  probability 
Assign  the  Z  value  corresponding  to  that  probability 
Standard  Normal  Inverse  Procedure 
Select  the  appropriate  range  of  standardized  Z  values 
Assign  the  probability  corresponding  to  that  range  of  Z  values 
Compute  "What-if'  Analysis  Procedure  (9  procedures,  each  for  a  given  number  of  sources) 
Declare  local  variables 

Initialize  variables  (including  proportions  (k,),  means,  and  variances  for  each  source) 
Determine  if  all  proportions  sum  to  1,  and  if  all  data  is  non-negative 
If  error  exists,  display  appropriate  message  and  exit  subroutine 
Use  a  loop  to  calculate  the  blend's  mean,  variance,  and  standard  deviation 
If  the  standard  deviation  is  zero,  then 
Set  the  probability  of  exceeding  maximum  safe  value  to  either  0  or  1 
Otherwise, 

Calculate  Z 

Call  the  Standard  Normal  Inverse  procedure  to  find  the  probability 
Optimize  Procedure  (9  procedures,  each  for  a  given  number  of  sources) 

Declare  local  variables 
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Initialize  absolute  bounds 
Initialize  source  means/  variances 

Initialize  Z,  current  conditions,  and  check  initial  centroid  feasibility 
Check  upper  bound  feasibility  (that  the  sum  of  upper  bounds  >  1) 

Set  initial  feasible  solution 

If  feasible,  then  initial  solution  =  centriod  point  (all  sources  set  to  1/NumSources) 
If  not,  then  set  all  sources  at  upper  bound  until  sum  =  1 

Compute  mean,  standard  deviation,  and  objective  function  value  of  initial  solution 
For  each  iteration: 

For  each  sample: 

Set  blend  proportions  (except  last  source)  randomly  within  current  boxmds 
Set  final  blend  proportion  as  (1-sum  of  others) 

Calculate  objective  function  value,  Mp 

Compare  Mp  to  Mp*,  replacing  the  optimal  if  Mp  is  less 

Repeat  for  all  samples 

After  all  samples  taken  (i.e.  end  of  iteration)  reduce  range  of  each  blending 
source  appropriately 
Repeat  for  all  iterations 

Calculate  final  values  to  send  as  output  (blend  mean,  standard  deviation,  objective 
function  value) 
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APPENDIX  10 

Code  for  Optimization  of  2  Sources 


Public  Sub  OptTwo (OptMean,  OptMax,  OptStandDev,  OptProb,  NumSources,  OptProport ( 

) 

As  Single) 

I  ***  this  procedure  is  the  code  for  optimizing  2  SOURCES  *** 


'  DECLARE  LOCAL  VARIABLES 
Dim  i  As  Integer 
Dim  i  As  Integer 
Dim  k  As  Integer 
Dim  Z  As  Single 

Dim  AbsLow(l  To  2)  As  Single 
Dim  AbsUpd  To  2)  As  Single 
Dim  Mud  To  2)  As  Single 
Dim  Yar  d  To  2)  As  Single 
Dim  OptimalMp  As  Single 
Dim  UBSum  As  Single 

Dim  BlendMean  As  Single 

Dim  BlendVar  As  Single 

Dim  CurrLow(l  To  2)  As  Single 

Dim  CurrRange(l  To  2)  As  Single 

Dim  CurrUpd  To  2)  As  Single 

Dim  CurrSoln(l  To  2)  As  Single 

Dim  CurrSolnMp  As  Single 

Dim  Sum  As  Single 

Dim  Feasible  As  Boolean 


index  for  sample  size 
index  for  iterations 
index  for  each  source 

Z  alpha  value  for  computing  max  value 

absolute  lower  bounds  for  each  source 
absolute  upper  bounds  for  each  source 
mean  output  for  each  source 
variance  for  each  source 
current  optimal  max  value 
error-checks  upper  bounds 

current  blend's  mean 
current  blend's  variance 
current  lower  bound  for  each  source 
current  feasible  range  for  each  source 
current  upper  bound  for  each  source 
current  solution's  proportions 
current  solution's  max  value 
determines  final  source ' s  proportion 
determines  initial  centroid  feasibility 


'  INITIALIZE  VARIABLES 
'  initialize  absolute  bounds 
For  k  =  1  To  NumSources 
AbsLow(k)  =  0 
N0X t  k 

AbsUpd)  =  frmSourceOne! txtSourceOneUB.Text 
AbsUp ( 2 )  =  f rmSourceTwo ! txtSour ceTwoUB . Text 
'  initialize  source  means/variances 

Mu(l)  =  frmSourceOne! txtSourceOneMean. Text 
Mu (2)  =  f rmSourceTwo ! txtSourceTwoMean. Text 

If  (frmSourceOne! optSourceOneSD. Value  =  True)  Then 
Yard)  =  (frmSourceOne! txtSourceOneVar. Text)  ^  2 
3-  s  0 

Var(l)  =  (frmSourceOne! txtSourceOneVar .Text) 

End  If 

If  (f rmSourceTwo !optSourceTwoSD. Value  =  True)  Then 
Var(2)  =  ( f rmSourceTwo !txtSourceTwoVar. Text)  2 
^■1  SC 

Var(2)  =  (frmSourceTwo! txtSourceTwoVar.Text) 

End  If  ,  . , 

'  initialize  Z,  current  conditions;  check  centroid  feasbility 
Feasible  =  True 
Z  =  ZAlpha (OptProb) 

UBSum  =  0 

For  k  =  1  To  NiomSources 
CurrLow(k)  =  AbsLow(k) 

CurrUp(k)  =  AbsUp (k) 

UBSum  =  UBSum  +  AbsUp (k) 

If  (AbsUp(k)  <  (1  /  NumSources))  Then 
Feasible  =  False 
End  If 
Next  k 

'  check  upper  bound  feasibility 

If  (UBSum  <  1)  Then  ^  „  .o 

MsgBox  "Upper  bounds  are  infeasible;  must  sum  to  at  least  1.  ,  4«, 

er  Bounds" 

Exit  Sub 


"Upp 


set  initial  feasible  solution;  compute  corresponding  obiective  function  val 


Sum  =  0 

BlendMean  =  0 

BlendVar  =  0 

For  k  =  1  To  NumSources 

If  (Feasible  =  True)  Then 

OptProport (k)  =  1  /  NumSources 
£  3.  s  6 

OptProport (k)  =AbsUp(k) 

Sum  =  Sum  +  OptProport (k) 

If  (Sum  >  1)  Then 

Sum  =  Sum  -  OptProport (k) 

OptProport (k)  =  1  -  Sum 
Sum  =  Sum  +  OptProport (k) 

End  If 
End  If 

BlendMean  =  BlendMean  +  (Mu(k)  *  OptProport (k) ) 

BlendVar  =  BlendVar  +  (  (OptProport  (k)  ''  2)  *  Var(k)) 

Ngx t  k 

OptimalMp  =  BlendMean  +  Z  *  Sqr (BlendVar) 

CONDUCT  OPTIMIZATION 
For  i  =  1  To  Numiterations 
For  i  =  1  To  SampleSize 
set  blend  proportions  randomly 
Sum  =  0 

For  k  =  1  To  (NumSources  -  1) 

CurrUp(k)  =  Min (CurrUp (k) ,  1  -  Sum) 

CurrRanae(k)  =  CurrUp (k)  -  CurrLow(k) 

CurrSoln(k)  =  CurrLow(k)  +  (CurrRanqe (k)  *  Rnd) 

Sum  =  Sum  +  CurrSoln(k) 

Next  k 

set  final  source  proportion  as  1  -  (Siam  of  others) 

CurrSoln (NumSources)  =  1  -  Sum 
calculate  obiective  function  value  (Mp) 

BlendMean  =  0 

BlendVar  =  0 

For  k  =  1  To  NumSources 

BlendMean  =  BlendMean  +  (Mu(k)  *  CurrSoln(k)) 

BlendVar  =  BlendVar  +  ( (CurrSoln (k)  ^  2)  *  Var(k)) 

Ngx t  k 

CurrSolnMp  =  BlendMean  +  Z  *  Sqr (BlendVar) 
compare  to  current  optimal 

If  CurrSolnMp  <=  OptimalMp  Then 
OptimalMp  =  CurrSolnMp 
For  k  =  1  To  NumSources 

OptProport (k)  =  CurrSoln (k) 

Next  k 
End  If 

Next  i  -  - 

after  all  samples  taken  (i.e.  end  of  iteration)  reduce  range  of  each 
blending  source  appropriately 
For  k  =  1  To  NumSources 

CurrLow(k)  =  Max  (OptProport  (k)  -  2  ''  (-1  -  1),  AbsLow(k)) 

CurrUp(k)  =  Min (OptProport (k)  +  2  ^  (-1  -  D /  AbsUp(k)) 

Next  k 
Next  1 

CALCULATE  VALUES  TO  SEND  AS  OUTPUT 
OptMean  =  0 
BlendVar  =  0 
For  k  =  1  To  NumSources 

If  OptProport (k)  <  0  Then 
OptProport (k)  =  0 
End  If 


OptMean  =  OptMean  +  (Mu(k)  *  OptProport (k) ) 
BlendVar  =  BlendVar  +  (  (OptProport  (k)  ''  2)  * 
Next  k 

OptStandDev  =  Sqr (BlendVar) 

OptMax  =  OptMean  +  Z  *  OptStandDev 

End  Sub 


Var (k) ) 
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Code  for  "What-if"  Analysis  for  2  Sources 


Public  Sub  WIComputeTwo (NumSources  As  Integer,  MaxSafeVal  As  Single, 
BMean  As  Single,  BStandDev  As  Single,  WIProb  As  Single) 

r  ***  THIS  PROCEDURE  IS  THE  CODE  FOR  CALCULATING 
'  A  "WHAT-IF"  ANALYSIS  FOR  2  SOURCES  *** 


'  Declare  Local  Variables 
Dim  i  As  Integer 
Dim  Ki(l  To  2)  As  Single 
Dim  Mu(l  To  2)  As  Single 
Dim  Var(l  To  2)  As  Single 
Dim  Z  As  Single 
Dim  NumericErr  As  Boolean 
Dim  ProportErr  As  Boolean 
Dim  MeanErr  As  Boolean 
Dim  VarErr  As  Boolean 
Dim  ProportSum  As  Single 
Dim  BVariance  As  Single 

'  Initialize  Variables 
ProportErr  =  False 
MeanErr  =  False 
VarErr  =  False 
ProportSum  =  0 

Ki(l)  =  frmSourceOne! txtSourceOneKi.Text 
Ki(2)  =  frmSourceTwo! txtSourceTwoKi.Text 

Mu(l)  =  frmSourceOne ! txtSourceOneMean. Text 
Mu (2)  =  frmSourceTwo ! txtSourceTwoMean. Text 

If  (frmSourceOne ! opt SourceOneSD. Value  =  True)  Then 
Var  ( 1 )  =  ( frmSourceOne !  txt Sour ceOneVar .  Text )  2 

Else 

Var(l)  =  (frmSourceOne! txtSourceOneVar. Text) 

End  If 

If  (frmSourceTwo! opt SourceTwoSD. Value  =  True)  Then 
Var (2)  =  (frmSourceTwo! txtSourceTwoVar. Text)  2 
E 1 S  S 

Var (2)  =  (frmSourceTwo! txtSourceTwoVar. Text) 

End  If 

'  Error-Check  Variables  (Probs  sum  to  1;  non-negative  values) 

For  i  =  1  To  2 

If  Ki(i)  <  0  Then 
ProportErr  =  True 
End  If 

If  Mu(i)  <  0  Then 
MeanErr  =  True 
End  If 

If  Var(i)  <  0  Then 
VarErr  =  True 
End  If 

ProportSum  =  ProportSum  +  Ki(i) 

Next  i 

'  If  Error  Exists,  Display  Appropriate  Message  and  Exit  Subroutine 
If  ( (ProportErr)  Or  (ProportSum  <>  1) )  Then  _  ^ 

MsgBox  "Proportions  must  be  non-negative  and  s\am  to  1.  Try  again.",  48,  'C 
annot  Compute" 

Exit  Sub 
End  If 

If  MeanErr  =  True  Then 

MsgBox  "Source  means  must  be  non-negative.  Try  again.",  48,  "Cannot  Comput 


'  Blend  Proportions 
'  Source  Mean 
'  Source  Variation 
'  Determines  Probability 
'  Detects  errors  in  Source  Data 
'  Detects  error  in  Blend  Proportion 
'  Detects  error  in  Source  Mean 
'  Detects  error  in  Source  Variation 
'  Sums  the  Blend  Proportions 
'  Blend  Variance 


e 


Exit  Sub 


i.lL«  «  0-1-0 


End  If 

If  VarErr ' =  True  Then 

MsqBox  "Source  variations  must  be  non-neqative .  Try  aqain.",  48, 
ompute" 

Exit  Sub 
End  If 

'  Calculate  Values 
For  i  =  1  To  2 

BMean  =  BMean  +  Ki(i)  *  Mu(i) 

BVariance  =  BVariance  +  (Ki(i)  2)  *  (Var(i)) 

Next  i 

BStandDev  =  Sqr (BVariance) 

If  BStandDev  =  0  Then 

If  BMean  >  MaxSafeVal  Then 
WIProb  =  1 
Else 

WIProb  =  0 
End  If 
Else 

Z  =  (MaxSafeVal  -  BMean)  /  BStandDev 
Call  SNormInv(Z,  WIProb) 

End  If 

End  Sub 


"Cannot  C 
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Validation  Data 


Validation  Data 


n  =  3 


Hi  =  50 
H2  =  40 
M3  =  45 


Oi^  =  4 
02^  =  36 
03^  =  12 


bi  =  1.00 

b2  =  1.00 
b3  =  1.00 


Mp*  =  52.16074 


experimental  runs 

1)  52.16074 

2)  52.16074 

3)  52.16074 

4)  52.16074 

5)  52.16074 


%  error  =  (52.16074  -  52.16074)  /  52.16074 
=  0.00% 


n  =  7 


Ml  =  50 
M2  =  40 
M3  =  45 
M4  =  51.5 
M5  =  47 
M6  =  51 
M7  =  42 


01^  =  4 
02^  =  36 
03^  =  12 
04^  =  0.01 
05^  =  20 
06^  =  1 
07^  =  27 


fci  =  1.00 
&2  =  1.00 
=  0.09 
&4  =  1.00 
bs  =  1.00 
fee  =  1.00 
b7=1.00 


Mp*  =  51.17793 


experimental  runs 

1)  51.19986 

2)  51.28685 

3)  51.23554 

4)  51.31067 

5)  51.28316 


%  error  =  (51.263216  -  51.17793)  /  51.17793 
=  0.17% 
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n  =  10 


pi  =  50 

=  4 

II 

o 

az^  =  36 

(J3  =  45 

032  =  12 

|J4  =  43 

042  =  17 

Ms  =  47 

052  =  21 

Me  =  41 

Ge^  =  14 

M7  =  42 

072  =  33 

Ms  =  49 

082  =  7 

M9  =  48 

092  =  11 

Mio  =  46 

Olo2  =  6 

experimental  runs 

1)  48.92244 

2)  49.30835 

3)  49.02160 

4)  48.92251 

5)  48.93250 


bi  =  1.00 

b2  =  1.00 
&3  =  1.00 
bi  =  1.00 
&5  =  1.00 
be  =  0.03 
&7=1.00 
bs  =  1.00 
b9  =  1.00 
bw  ~  1.00 

Mp*  =  48.92204 


%  error  =  (49.02148  -  48.92204)  /  48.92204 
=  0.20% 
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